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Abstract 

In this paper we describe the basic features of the Bergm package for the 
open-source R software which provides a comprehensive framework for Bayesian 
analysis for exponential random graph models: tools for parameter estimation, 
model selection and goodness-of-fit diagnostics. We illustrate the capabilities of 
this package describing the algorithms that drive the package through a tutorial 
analysis of two well-known network datasets. 

1 Introduction 

The R package Bergm implements Bayesian analysis for Exponential Random Graph 
Models (ERGMs) (Wasserman and Pattison (1996); Robins et al. (2007)) using the 
methods described by Caimo and Friel (2011) and Caimo and Friel (2012). The package 
provides a comprehensive framework for Bayesian inference and model selection using 
Markov chain Monte Carlo (MCMC) algorithms. It can also supply graphical Bayesian 
goodness-of-fit procedures that address the issue of model adequacy. Although com- 
putationally intensive, the package is simple to use and represents an attractive way 
of analyzing network data as it offers the advantange of a complete probabilistic treat- 
ment of uncertainty. Bergm is based on the ergm package developed by Hunter et al. 
(2008b) and therefore it makes use of the same model set-up and network simulation 
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algorithms. The Bergm package has been continually improved in terms of speed per- 
formance over the last two years and one of the purposes of this paper is to highhght 
these improvements. We feel that this package now offers the end-user a feasible option 
for carrying out Bayesian inference for exponential random graphs. 

Two well-known network datasets will be used throughout this tutorial for illus- 
trative purposes: the first is the Kapferer Tailor Shop dataset (Kapferer, 1972) whose 
directed edges represent work interactions in a tailor shop in Zambia (then North- 
ern Rhodesia) and nodai attributes refer to the job status. The second network is 
Zachary's karate club (Zachary, 1977) which represents the undirected social network 
graph of friendships between 34 members of a karate club at a US university in the 
1970s. Figure 1 displays the graphs of these two networks. 

In this paper we describe how to instali and load Bergm (Section 2) providing a brief 
summary of what Bayesian ERGMs are (Section 3). Sections 4, 5, and 6 overview the 
algorithms and the functions used to produce posterior estimates for the parameters, 
Bayesian goodness-of-fit procedures and model selection respectively. This paper does 
not provide an exhaustive description of all the functionality and options available, and 
more Information about the commands and methods mentioned are available through 
the R help system within the package. 




Figure 1: Kapferer Tailor Shop directed graph (right) and Zachary's karate club undi- 
rected graph (right). 
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2 Getting Bergm 



The Bergm package can be obtained and loaded in R using the following commands: 

> instali .packagesC "Bergm") 

> libraryC "Bergm") 

Since Bergm depends on ergm (which in turn depends on network), coda, and mvt- 
norm, installing the package will automatically load all the dependencies. Ali of 
these packages are available on the Comprehensive R Archive Network (CRAN) at 
http://CRAN.R-project.org. 

The results presented in this paper have been obtained using R version 2.13.2 on 
a Mac using Bergm version 2.1; ergm version 2.4-3; network version 1.6; coda version 
0.14-6; and mvtnorm version 0.9-999. 



3 Bayesian exponential random graphs 

The Bayesian approach to statistical problems is probabilistic. Inference is based on the 
posterior distribution which is the conditional probability of the unknown quantities 
given the observed ones. The posterior extracts the information in the data and provide 
a complete summary of the uncertainty about the unknowns. 

In the ERGM context (see Wasserman and Pattison (1996) and Robins et al. 
(2007)), the purpose of Bayesian inference is to learn about the posterior distribu- 
tion of the model parameters 6 of an observed graph y n nodes: 

where s{y) is a known vector of sufficient network statistics (Figure 2) (Morris et al., 
2008), p{0) is a prior distribution placed on 6, z{0) is the likelihood normalizing con- 
stant, and p{y) is the model evidence. Equation (1) provides a probabilistic statement 
about how likely parameter values are after observing the data y. The likelihood 
p{y\0) is translated into a proper probability distribution that can be summarised by 
computing expected values, Standard deviations, quantiles, etc. 

Unfortunately the posterior distribution (1) is doubly-intractable as both z{6) and 
p{y) cannot be evaluated analytically (Koskinen, 2004). This makes the use of Standard 
MCMC procedures infeasible. 
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Figurc 2: Some of the most used configurations for directed (a) and undirected (b) 
graphs. 
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In order to carry out Bayesian inference for ERGMs, the Bergm package makes use 
of a combination of Bayesian algorithms and MCMC techniques. The exchange algo- 
rithm circumvents the problem of computing the normahzing constants of the ERGM 
hkehhoods, while the use of muhiple chains interacting with each others (population 
MCMC approach) by means of adaptive direction samphng is able to speed up the 
computations and improve chain mixing quite significantly. 

4 Bayesian parameter estimation 

In order to estimate (1), the Bergm package uses the exchange algorithm described in 
Section 4.1 of Caimo and Friel (2011) to sample from the foUowing distribution: 

p{0\y\0\y)^p{y\e)p{e)e{e'\e)p{y'\e') (2) 

where p{'ij\6') is the hkehhood on which the simulated data y' are defined, e{6'\0) is 
any arbitrary proposal distribution for the augmented variable 6' . As we will see in 
the next section, this proposal distribution is set to be a normal centered at 6. 

At each MCMC iteration, the exchange algorithm consists of a Gibbs update of 
6' followed by a Gibbs update of which is drawn from the p{-\6') via an MCMC 
algorithm (Hunter et al., 2008b). Then a deterministic exchange or swap from the 
current state 6 to the proposed new parameter 6' . This deterministic proposal is 
accepted with probability: 

. / qg{,/)p{e')h{e\e')qe,{y) z{e)z{e') \ 
" v ' qo{y)pmmO)qe'{y') z{0)z{0') ) ' 

where qg and g^/ indicates the unnormalised likelihoods with parameter 6 and O', 
respectively. Notice that all the normalising constants cancel above and below in the 
fraction above, in this way avoiding the need to calculate the intractable normalising 
constant. 

The exchange algorithm is implemented by the bergm function in the following 
way: 



for i = 1, . . . , 

1. generate 6' from e{-\d) 
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2. simulate iJ fromp(-|0') 



3. update ^ 0' 



with probability 



log(a) = min O, [0 - 6»']* [s{y') - s{y)] + log 



\p{e') 



) 



(3) 



V{0) 



end for 



where s[y) is the observed vector of network statistics and s(y') is the simulated vector 
of network statistics. 

4.1 Block-update sampler 

Step 1 of the algorithm consists in generating 6' from some proposal distribution 
within each iteration. Bergm uses a block-update sampler with normal proposal to 
simultaneously update of the parameter values in the MCMC chain: 



Typically, tuning the proposal distribution e from which O' is drawn represents the 
crucial part of the algorithm since a poor tuning of the proposal parameter Sg can 
slow down the chain's mixing rate and therefore the algorithm can take a very long 
time to converge to the stationary posterior density. 

4.2 Parallel adaptive direction sampler 

In order to improve mixing a parallel adaptive direction sampler (ADS) is considered: 
at each i-th iteration of the algorithm we have a coUection of H different chains in- 
teracting with one another. By construction, the state space consists of {^i, . . . , Oh} 
with target distribution p{Oi\y) ® ■ ■ ■ ®p{OH\y)- A parallel ADS move consists of gen- 
erating a new value 0'^^ from the difference of two parameters 0^^ and 0^^ (randomly 
selected from other chains) multiplied by a scalar term 7 which is called parallel ADS 
move factor plus a random term e called parallel ADS move parameter (Figure 3) 
which is equivalent to the block-update sampler defined in (4). The algorithm can be 
summarised as follows: 



e{0'\0)^M{0,i:,). 



(4) 
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for i = 1, . . . ,iV 

for h = l,...,H 

1 . generate hi and /i2 such that hi ^ h2 ^ h 

2. generate 0^ from 'y{6hi - Ohi) + 

3. simulate i/ from p{-\6'f^) 

4. update 6h ^ 0'^ with probability 



min O, [eh-e'^X[s{y')-s{y)]+\og 



end for 
end for 



(5) 




Figure 3: The parallel ADS move of the current state (darker blue dot) consists of 
generating a new parameter value along the direction made by the difference of two 
randomly sampled parameter states (hght blue dots) belonging to different chains plus 
a random term e. 
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4.3 Kapferer tailor shop network 

Consider the Kapferer Tailor Shop network and a 3-dimensional model including the 
following network statistics: edges (edges), mutual edges (mutual) and cyclic triples 
(ctriple) involving nodes with the same job status, where the job status is represented 
by a categorical variable of 8 levels: 

mutual Y.i<j VijVji 

ctripleC job") J2i<j<kyijyjkyki where i, j, k have the same job status 
The format of the model specification is the same of an ergm formula: 

> formula <- y ~ edges + mutual + ctripleC job") 

Then we can use the bergm function to sample from the posterior distribution using the 
MCMC algorithm described above. In this example we use the parallel ADS procedure. 
By default, the number of chains is set as twice the number of dimensions of the model. 
It is possible to choose a different number of chains by using the argument nchains. 
In order to perform the block-site update described in Section 4.1 it is necessary to 
set nchains = 1. For each chain, we can then set the number of burn-in iterations 
(burn. in) and the number of iterations after the burn-in (main. iters). The number 
of iterations used to simulate a network y' at each iteration is defined by the argument 
aux. iters. 

> post.est <- bergm(f ormula, 



+ burn.in=200, 

+ gamma=0 . 7 , 

+ main. iters=1000, 

+ aux. iters=25000) 



The population MCMC with parallel ADS move is the default procedure of the bergm 
function. The total number of iterations (eg the size of the posterior sample) is nchains 
X main. iters. The proposal covariance structure is defined by the argument 
sigma, epsilon which is set to be a diagonal matrix with every diagonal entry equal to 
a small number. In many cases, good mixing of the chain is ensured by a sensible tuning 
of the parallel ADS move factor gamma and therefore the argument sigma, epsilon can 
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be generally left at its default setting. As said above, parallel ADS is adopted as the 
default procedure but it is automatically disabled in the case of uni-dimensional models 
where the block-update sampler is used and the argument gamma is used to tune the 
varian ce of the normal proposal distribution e. 

After completing the estimation, post.est is an object of the class bergm and 
contains a list of attributes among which are the real and CPU time (in seconds) 
taken by the estimation process: 

> post . est$time 

user system elapsed 
99.083 0.759 100.659 

It is possible to visualise the results of the MCMC estimation by using the bergm . output 
function which is based on the coda package (Plummer et al., 2006) which is an R pack- 
age for MCMC output analysis and diagnostics. 

> bergm . output (post . est , lag . max=50) 

MCMC results for Model: y ~ edges + mutual + ctripleC job") 
Posterior mean: 





thetal 


(edges) 


theta2 (mutual) 


thetaS 


(ctriple .job) 


Chain 1 


-3, 


.4256853 


3 . 7647452 




0.8168465 


Chain 2 


-3, 


.4091174 


3 . 7447376 




0.8567871 


Chain 3 


-3, 


.4140165 


3.7377072 




0.8532598 


Chain 4 


-3, 


.4366544 


3.7868822 




0.8585445 


Chain 5 


-3, 


.4201932 


3.8034972 




. 8340857 


Chain 6 


-3, 


.4009597 


3.7021253 




0.8433665 


Posterior 


• sd: 












thetal 


(edges) 


theta2 (mutual) 


thetaS 


(ctriple .job) 


Chain 1 


0, 


. 2046574 


0.4287769 




0.1882377 


Chain 2 


0, 


.2291025 


0.4960458 




0.1829550 


Chain 3 


0, 


. 1903504 


0.4453100 




0.1700874 


Chain 4 


0, 


.2673414 


0.5266250 




0.1655076 
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Chain 5 0.2474045 0.4652013 0.1543426 

Chain 6 0.2696806 0.4775102 0.1949814 







Acceptance rate: 


Chain 


1 


0.213 


Chain 


2 


0.267 


Chain 


3 


0.246 


Chain 


4 


0.224 


Chain 


5 


0.260 


Chain 


6 


0.237 



Overall posterior density estimate: 

thetal (edges) theta2 (mutual) theta3 (ctriple . job) 

Post. mean -3.4177711 3.7499491 0.8438150 

Post. sd 0.2373836 0.4768709 0.1782518 

Overall acceptance rate: 0.241166666666667 

The output above shows the results of the MCMC estimation: posterior means and 
Standard deviations, and acceptance rates for every chain in the population and for 
the overall chain. Notice that the posterior summaries of each chain are consistent 
with each other. Figure 4 displays the MCMC diagnostic plots produced by the 
bergm . output function. In this example, we notice a low density effect expressed 
by the negative value of the posterior mean of the density effect parameter (edges) 
combined with the positive mutuality and transitivity within people having the same 
job status expressed by the mutual edge and cyclic triple parameters respectively. The 
overall acceptance rate is around 24% and the autocorrelation is negligible after lag 
50. 
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MCMC output for Model: y ~ edges + mutual + ctriple("job") 



I 1 1 1 

-4.0 -3.5 -3.0 -2.5 



CO 



o 




1 — I — I — \ — I — I — r 

o 2000 4000 6000 



LO 

d 



o 
d 



1 — I — I — I — I — r 
O 10 20 30 40 50 



S-i (edges) 



Iterations 



Lag 




o 
iri 



o 



o 
c6 



o 




1 — I — I — \ — I — I — r 

o 2000 4000 6000 



o 
d 




1 — I — I — I — I — r 

o 10 20 30 40 50 



82 (mutual) 



Iterations 



Lag 




I 1 1 

0.5 1.0 1.5 



CO 

d 



CM 

d 




1 1 1 1 1 1— 

O 2000 4000 6000 



o 
d 




T 1 1 1 1 T 

O 10 20 30 40 50 



83 (ctriple.job) 



Iterations 



Lag 



Figure 4: MCMC diagnostics for the overall chain. 
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5 Bayesian goodness-of-fit diagnostics 



The bgof function provides a useful tool for assessing Bayesian goodness-of-fit so as to 
examine the fit of the data to the posterior model obtained by the bergm function. The 
observed network data y is compared with a set of networks Ui, ■ ■ ■ , Us simulated 
from S independent reahsations 0i, 62, . . . ,0s of the posterior density estimate. This 
comparison is made in terms of high- level characteristics g{-) such as higher degree 
distributions, etc. (see Hunter et al. (2008a)). The algorithm can be summarised as 
foUows: 



for i = 1, . . . , S 

1. sample 6i from the estimate of p{0\y) 

3. simulate f rom j9(-|0i) 

4. calculate g{yi) 

end for 



For example, the code below is used to compare the Kapferer Tailor Shop network 
with a series of networks simulated from 100 random reahsations (sample . size) of 
the estimated posterior distribution post . est using 50, 000 iterations (aux. iters) for 
the network simulation step. The bgof function may take a few seconds to run and, 
at the end of the execution, it will automatically plot the results as shown in Figure 5. 

> bgof (post . est , 

+ sample . size=100, 

+ aux. iters=50000, 

+ directed=TRUE, 

+ n.ideg=20, 

+ n.odeg=20, 

+ n.dist=10, 

+ n.esp=15) 

The set of statistics used for the comparison of directed networks includes the in-degree 
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distribution, the out-degree distribution, the minimum geodesic distance distribution 
and the edgewise shared partner distribution. The arguments n . ideg, n . odeg, n . dist, 
and n.esp indicates the number of boxplots to plot for each distribution respectively. 




123456789 02468 10 13 

minimum geodesic distance edge-wise sliared partners 



Figure 5: Bayesian goodness-of-fit diagnostics. The red hne displays the goodness of 
fit statistics for the observed data together with boxplots of goodness of fit statistics 
based on 100 simulated networks from the posterior distribution. 

In Figure 5 we see, based on the various goodness of fit statistics, that the networks 
simulated from the posterior distribution are in reasonable agreement with the observed 
network. We can therefore conclude that the data is a reasonable fit to the model, 
despite its simplicity. 
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6 Bayesian model selection 



An important problem in statistical analysis is the choice of an optimal model from a 
set of a priori competing models. In the ERGM context, tliis task translates into the 
choice of which subset of network statistics should be included into the model. 

Let m indicate a particular model from a set of competing models with correspond- 
ing parameters 0. Following the Bayesian paradigm, interest focuses on exploring the 
posterior distribution, 

plrrijdly) (X p{y\m,0) plOlm) p{m) (6) 

where p{6\m) and p(m) are prior distributions within model m, and on model m, 
respectively. The reversible Jump Markov chain Monte Carlo (RJMCMC) algorithm 
(Green, 1995) was designed to explore this type of posterior distribution across the 
joint model and parameter space. It is therefore a type of MCMC algorithm that 
allows one to jointly explore the uncertain between and within models. 

This approach is very appealing since it relies exclusively on probabilistic consid- 
erations but is very challenging from a computational viewpoint. As stated above, 
the intractability of the likelihood normalising constant z{0) in (1) renders Standard 
RJMCMC techniques infeasible. However, the exchange algorithm used for parameter 
estimation can be easily generalised so as to include model indicators. 

The auto-RJ exchange algorithm described in Caimo and Friel (2012) represents a 
trans-dimensional RJMCMC extension of the exchange algorithm involving an inde- 
pendence sampler based on a distribution fitting a parametric density approximation 
to the within-model posterior. This approach overcomes the issue of the likelihood 
intractability sampling from: 

p{0', 0, m' , m, l/lu) (X p{y\0, m)p{0\m)p{m)w{0' \m')h{m' \m)p{'i/ \0' , m') (7) 

where m and m' are two competing models, p{y\0,m) and p{'i/\0' ,m') are the two 
likelihoods for the data y under model m and the simulated data y' under model 
m' respectively, p{0\m) and p{m) are the priors for the parameter and the respective 
model m, w{0'\m') is a within-model proposal (independence sampler) which fit a 
parametric density approximation to the model posteriors and h{m'\m) is a between- 
model proposal. Notice that the marginal distribution for for 0' and m' in (7) is the 
target distribution of interest p{0,m\y). 
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The bergmS function implements the auto-RJ exchange algorithm which consists 
of two parts: an offline step and an online step. In the first step, samples from the 
posterior p{d\y, m) are gathered from each competing model using the bergm function 
and then approximated by normal distributions M{iimi S^) determined by the first 
and second moments from a sample from the model. 

The second step (online run) of the algorithm consists of a Gibbs update of m' 
followed by a Gibbs update of 6' which is generated via the independence sampler 
w{6'\m') ~ Af{fim', Sm'). This is followed by a Gibbs update of y' which is generated 
from p{-\0',m'). Then a deterministic exchange move from a current state {0,m) to 
the proposed new state {d', m') is accepted with probability: 

. J. qe,m{'}/) q9',m'{y) p{0'\m')p{m') w{e\fim,'Sr 

a = mm < 1, — - — -— — r — - — ^ 

[ leMV) qe',m'[y!) p[0\m) p{m) w{e'\fi„,> ,^„,>) 

where g^^m and qe',m' indicates the unnormalised likelihoods under model m with pa- 
rameter 6 and under model m' with parameter 6' respectively. 

The structure of the bergmS function can be described in the following way: 



for i = 1, . . . , 

1. generate m' from /i(-|m) 

2 . generate O' ~ J^^fim', S^') 

3. simulate i/ f rom p(-|0', m') 

4. update {0,m) ^ {O', m') with probability: 

log(a) = min (o, 0\sm{lJ) - Sm{y)] + 0\sm'{lJ) - Sm'(z/)] + log 
end for 



p{0') w{0\flra,'Sr, 

p{0) w{0'\ii^,,% 



where Sm{y) and Sm'{y) are the observed vectors of network statistics under model 
m and m' respectively, and Sm{lJ) and Sm'{'}/) are the simulated vector of network 
statistics under model m and m' respectively. The reader is referred to Caimo and 
Friel (2012) for much details on this algorithm. 
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6.1 Karate club network 



Consider the Karate club network, we propose three models to fit the data using a set of 
new specification statistics introduced by Snijders et al. (2006): geometrically weighted 
edgewise shared partners (gwesp) and geometrically weighted degrees (gwdegree): 



where the scale parameters 0^ = 0.2 and 0„ = 0.8. Dk{y) is the number of pairs that 
have exactly k common neighbours and EPk{y) is the number of connected pairs with 
exactly k common neighbours. The specification of these models requires the creation 
of a list of formulas: 

> formulae <- c(y ~ edges + gwesp(0.2,f ixed=TRUE) , 

+ y ~ edges + gwdegree(0.8,f ixed=TRUE) , 

+ y ~ edges + gwesp(0.2,f ixed=TRUE) 

+ + gwdegree(0.8,f ixed=TRUE)) 

The bergmS command is then used to carry out the algorithm. To do this we have to 
specify several arguments for both the offline and online step. 

The offline run consists of running the bergm function for each of the models pro- 
posed. Therefore we set some arguments main.iters, burn.ins, gammas which are 
vectors containing values for bergm arguments: main. iters, burn. in, gamma for each 
competing model. 

The argument iters refers to the number of iterations used for the online run. 
The number of MCMC steps used for network simulation is specified as usual by the 
argument aux. iters and this will be used in both the offline and the online step. The 
command below should take around 10 minutes depending on the CPU speed of the 
Computer. 

> mod.sel <- bergmS (formulae, 

+ iters=25000, 

+ aux. iters=15000, 

+ main. iters=rep (700,3) , 

+ burn. ins=rep( 100,3) , 

+ gammas=c(l,l,0.8)) 



gwdegree 



gwesp 
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The bergmS . output function produces the MCMC diagnostics for each competing 
model explored by the MCMC algorithm. Figure 6 and 7 displays the plots regarding 
the posterior model and parameter density estimate respectively. 

> best.mod <- bergmS. output (mod. sel) 



BEST MODEL 



Model 1: y ~ edges + gwesp(0.2, fixed = TRUE) 



Posterior parameter est 

thetal (edges) 

theta2 (gwesp . f ixed . O . 2) 

Within-model acceptance 



mate : 

Post . mean: Post . sd: 

-3.2574625 0.3278196 

1.1008261 0.2515162 

rate: 0.26 



Model 3: y ~ edges + gwesp (0.2, fixed = TRUE) + gwdegree (O . 8 , fixed = TRUE) 



Posterior parameter est 

thetal (edges) 

theta2 (gwesp.f ixed.0.2) 

theta3 (gwdegree) 



mate : 

Post. mean: Post. sd: 

-3.4916584 0.5146103 

1.1824831 0.2737858 

0.4783638 0.5852124 



Within-model acceptance rate: 0.14 



BF_13 = 13.4508670520231 



Between-model acceptance rate: 0.03 

In the results above we have the posterior parameter estimates for two of the compet- 
ing models (Model 2 has not been visited through the MCMC runs) with respective 



17 



within-model acceptance rates and an estimate of the Bayes Factor (about 13) for the 
comparison between Model 1 and Model 3 which makes clear that tliere is evidence 
that Model 1 is the best model of the set. 

After running the command bergmS . output, it is possible to perform a Bayesian 
goodness-of-fit tests. In this case, since the observed network is undirected, the set 
of high-level statistics include the degree distribution in place of the in-degree and 
out-degree distributions. 

> bgof (best .mod, 

+ aux. iters=30000, 

+ n.deg=20, 

+ n.dist=10, 

+ n.esp=15) 

In this example, the observed data appears to be a reasonable fit to the posterior 
distribution of the model selected (Model 1), based on the goodness of fit geodesic and 
the shared partners statistics displayed in Figure 8. 



7 Discussion 

The software package Bergm aims to help researchers and practitioners in two ways. 
Firstly, it is currently the only package for R that provides a simple and complete 
range of tools for conducting Bayesian analysis for exponential random graph models. 
Secondly, Bergm makes available a platform that can be easily customised, extended, 
and adapted to address different requirements. 

The software package is under continual development and it is far from finished. 
The main limitation of the software is its computational cost which makes it unsuitable 
for managing network graphs larger than hundreds of nodes, however it is perfectly 
suited for networks involving up to a hundred nodes. However an important improve- 
ment in terms of computational time and efficiency will be done by turning some of 
the R functions into C functions integrated with the ergm package. We expect that 
will yield further reductions in computional run time. 

Future versions of the Bergm package will address several issues including Bayesian 
analysis of Curved Exponential Random Graph Models (Hunter and Handcock, 2006) 
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Figure 6: MCMC diagnostics: posterior model probabilities. 
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Figure 7: MCMC diagnostics: posterior parameter probabilities. 
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Figure 8: Bayesian goodness-of-fit diagnostics. 
and exponential random graph models with missing data. 
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